{ *********************************************************** }
{ *                    ksTools Library                      * }
{ *       Copyright (c) Sergey Kasandrov 1997, 2009         * }
{ *       -----------------------------------------         * }
{ *         http://sergworks.wordpress.com/kstools          * }
{ *********************************************************** }

unit TestFFT;
{

  Delphi DUnit Test Case
  ----------------------
  This unit contains a skeleton test case class generated by the Test Case Wizard.
  Modify the generated code to correctly setup and call the methods from the unit
  being tested.

}

interface

uses
  TestFramework, Windows, Forms, Dialogs, Controls, Classes, SysUtils, Variants,
  Graphics, Messages, ksMath;

type
  PComplexArray = ^TComplexArray;
  TComplexArray = packed array[0..$FFFF] of TksComplex;
  TDynComplexArray = packed array of TksComplex;

  PRealArray = ^TRealArray;
  TRealArray = packed array[0..$FFFF] of Extended;
  TDynRealArray = packed array of Extended;

  TRealFunc = reference to function(N: Integer): Extended;

type
  TTestFFT = class(TTestCase)
  private
    FFunc0, FFunc1, FFunc2: TRealFunc;
    procedure TestComplexA(Data: PComplexArray; Count: Integer; Sign: Integer);
    procedure TestComplexB(Data: PComplexArray; Count: Integer);
    procedure TestRealA(Data: PRealArray; Count: Integer; Sign: Integer);
    procedure TestRealB(Data: PRealArray; Count: Integer; Sign: Integer);
    procedure TestRealC(Data: PRealArray; Count: Integer);
    procedure TestRealF(Func: TRealFunc; Count: Integer; Sign: Integer);
    procedure TestTwoA(Data1, Data2: PRealArray; Count: Integer);
    procedure TestTwoF(Func1, Func2: TRealFunc; Count: Integer);
    procedure TestTwoG(Func1, Func2: TRealFunc; Count: Integer);
    procedure TestSinF(Func: TRealFunc; Count: Integer);
    procedure TestCos1A(Data: PRealArray; Count: Integer);
    procedure TestCos1F(Func: TRealFunc; Count: Integer);
    procedure TestCos2A(Data: PRealArray; Count: Integer);
    procedure TestCos2F(Func: TRealFunc; Count: Integer);
    procedure TestCorrA(Data1, Data2: PRealArray; Count: Integer);
    procedure TestCorrF(Func1, Func2: TRealFunc; Count: Integer);
    procedure TestAutoCorrA(Data: PRealArray; Count: Integer);
    procedure TestAutoCorrF(Func: TRealFunc; Count: Integer);
    procedure TestSpectrumA(Data: PRealArray; Count: Integer);
    procedure TestSpectrumF(Func: TRealFunc; Count: Integer);
  protected
    procedure SetUp; override;
  published
    procedure TestComplexA0;
    procedure TestComplexA1;
    procedure TestComplexB0;
    procedure TestRealA0;
    procedure TestRealA1;
    procedure TestRealB0;
    procedure TestRealB1;
    procedure TestRealC0;
    procedure TestRealF0;
    procedure TestTwoA0;
    procedure TestTwoF0;
    procedure TestTwoG0;
    procedure TestCorrA0;
    procedure TestCorrF0;
    procedure TestAutoCorrA0;
    procedure TestAutoCorrF0;
    procedure TestSpectrumA0;
    procedure TestSpectrumF0;
    procedure TestSinF0;
    procedure TestCos1A0;
    procedure TestCos2A0;
  end;

implementation

// Discrete Fourier Transform of a Complex Function
function ComplexDFT(Data: PComplexArray; Count, Sign: Integer): TDynComplexArray;
var
  M, N: Integer;
  Arg: TksComplex;

begin
  SetLength(Result, Count);
  for M:= 0 to Count - 1 do begin
    Result[M]:= 0;
    for N:= 0 to Count - 1 do begin
      Arg.Re:= 0;
      Arg.Im:= (2 * Pi * N * M) / Count;
      if Sign < 0 then Arg.Im:= - Arg.Im;
      Result[M]:= Result[M] + Data^[N] * TksComplex.Exp(Arg);
    end;
  end;
end;

// Discrete Fourier Transform of a Real Function
function RealDFT(Data: PRealArray; Count, Sign: Integer): TDynComplexArray; overload;
var
  M, N: Integer;
  Arg: TksComplex;

begin
  SetLength(Result, Count);
  for M:= 0 to Count - 1 do begin
    Result[M]:= 0;
    for N:= 0 to Count - 1 do begin
      Arg.Re:= 0;
      Arg.Im:= (2 * Pi * N * M) / Count;
      if Sign < 0 then Arg.Im:= - Arg.Im;
      Result[M]:= Result[M] + Data^[N] * TksComplex.Exp(Arg);
    end;
  end;
end;

function RealDFT(Func: TRealFunc; Count, Sign: Integer): TDynComplexArray; overload;
var
  M, N: Integer;
  Arg: TksComplex;

begin
  SetLength(Result, Count);
  for M:= 0 to Count - 1 do begin
    Result[M]:= 0;
    for N:= 0 to Count - 1 do begin
      Arg.Re:= 0;
      Arg.Im:= (2 * Pi * N * M) / Count;
      if Sign < 0 then Arg.Im:= - Arg.Im;
      Result[M]:= Result[M] + Func(N) * TksComplex.Exp(Arg);
    end;
  end;
end;

function SinDFT(Func: TRealFunc; Count: Integer): TDynRealArray; overload;
var
  M, N: Integer;
  Arg: Extended;

begin
  SetLength(Result, Count);
  for M:= 0 to Count - 1 do begin
    Result[M]:= 0;
    for N:= 0 to Count - 1 do begin
      Arg:= (Pi * N * M) / Count;
      Result[M]:= Result[M] + Func(N) * Sin(Arg);
    end;
  end;
end;

{ === First Form of Cosine Transform === }

function Cos1DFT(Func: TRealFunc; Count: Integer): TDynRealArray; overload;
var
  K, N: Integer;
  Arg: Extended;

begin
  SetLength(Result, Count + 1);
  for K:= 0 to Count do begin
    if Odd(K) then
      Result[K]:= 0.5 * (Func(0) - Func(Count))
    else
      Result[K]:= 0.5 * (Func(0) + Func(Count));
    for N:= 1 to Count - 1 do begin
      Arg:= (Pi * N * K) / Count;
      Result[K]:= Result[K] + Func(N) * Cos(Arg);
    end;
  end;
end;

function Cos1DFT(Data: PRealArray; Count: Integer): TDynRealArray; overload;
begin
  Result:= Cos1DFT(
    function(N: Integer): Extended
    begin
      Result:= Data[N];
    end,
    Count);
end;

{ === Second Form of Cosine Transform === }

function Cos2DFT(Func: TRealFunc; Count, Sign: Integer): TDynRealArray; overload;
var
  K, N: Integer;
  Arg: Extended;

begin
  SetLength(Result, Count);
  for K:= 0 to Count - 1 do begin
    Result[K]:= 0;
    for N:= 0 to Count - 1 do begin
      Arg:= (Pi * (N + 0.5) * K) / Count;
      Result[K]:= Result[K] + Func(N) * Cos(Arg);
    end;
  end;
  if Sign < 0 then Result[0]:= 0.5 * Result[0];
end;

function Cos2DFT(Data: PRealArray; Count, Sign: Integer): TDynRealArray; overload;
begin
  Result:= Cos2DFT(
    function(N: Integer): Extended
    begin
      Result:= Data[N];
    end,
    Count,
    Sign);
end;

function ComplexArray(Data: PComplexArray; Count: Integer): TDynComplexArray;
begin
  SetLength(Result, Count);
  Move(Data^, Result[0], Count * SizeOf(TksComplex));
end;

function RealArray(Data: PRealArray; Count: Integer): TDynRealArray; overload;
begin
  SetLength(Result, Count);
  Move(Data^, Result[0], Count * SizeOf(Extended));
end;

function RealArray(Func: TRealFunc; Count: Integer): TDynRealArray; overload;
var
  M: Integer;

begin
  SetLength(Result, Count);
  for M:= 0 to Count - 1 do
    Result[M]:= Func(M);
end;

// circular correlation calculated directly
function DirectRealCorr(Func1, Func2: TRealFunc;
                        Count: Integer): TDynRealArray; overload;
var
  M, N, I: Integer;

begin
  SetLength(Result, Count);
  for M:= 0 to Count - 1 do begin
    Result[M]:= 0;
    for N:= 0 to Count - 1 do begin
      I:= N + M;
      if I >= Count then I:= I - Count;
      Result[M]:= Result[M] + Func1(I) * Func2(N);
    end;
  end;
end;

function DirectRealCorr(Data1, Data2: PRealArray;
                        Count: Integer): TDynRealArray; overload;
begin
  Result:= DirectRealCorr(
    function(N: Integer): Extended
    begin
      Result:= Data1[N];
    end,
    function(N: Integer): Extended
    begin
      Result:= Data2[N];
    end,
    Count
  );
end;

{
function DirectRealCorr(Data1, Data2: PRealArray;
                        Count: Integer): TDynRealArray; overload;
var
  M, N, I: Integer;

begin
  SetLength(Result, Count);
  for M:= 0 to Count - 1 do begin
    Result[M]:= 0;
    for N:= 0 to Count - 1 do begin
      I:= N + M;
      if I >= Count then I:= I - Count;
      Result[M]:= Result[M] + Data1[I] * Data2[N];
    end;
  end;
end;
}

// circular correlation calculated directly
function DirectRealAutoCorr(Func1: TRealFunc;
                            Count: Integer): TDynRealArray; overload;
var
  M, N, I: Integer;

begin
  SetLength(Result, Count);
  for M:= 0 to Count - 1 do begin
    Result[M]:= 0;
    for N:= 0 to Count - 1 do begin
      I:= N + M;
      if I >= Count then I:= I - Count;
      Result[M]:= Result[M] + Func1(I) * Func1(N);
    end;
  end;
end;

function DirectRealAutoCorr(Data1: PRealArray;
                            Count: Integer): TDynRealArray; overload;
begin
  Result:= DirectRealAutoCorr(
    function(N: Integer): Extended
    begin
      Result:= Data1[N];
    end,
    Count
  );
end;

const
  ComplexArr0: packed array[0..3] of TksComplex = (
   (Re:0.5; Im:1.5), (Re:2.2; Im:1.4), (Re:3.3; Im:2.4), (Re:1.2; Im:2.3)
  );
  ComplexArr1: packed array[0..3] of TksComplex = (
   (Re:0.5; Im:0.5), (Re:0.25; Im:0.75), (Re:0.75; Im:0.5), (Re:0.5; Im:0.25)
  );
  ComplexArr2: packed array[0..7] of TksComplex = (
   (Re:0.5; Im:1.5), (Re:2.2; Im:1.4), (Re:3.3; Im:2.4), (Re:1.2; Im:2.3),
   (Re:0.5; Im:0.5), (Re:0.25; Im:0.75), (Re:0.75; Im:0.5), (Re:0.5; Im:0.25)
  );

const
  RealArr0: packed array[0..3] of Extended =
    (0.5, 1.1, 2.2, 3.3);
  RealArr1: packed array[0..7] of Extended =
    (0.8, 1.2, 0.2, 2.5, 1.4, 1.9, 2.6, 3.1);
  RealArr2: packed array[0..7] of Extended =
    (0.2, 3.4, 1.9, 2.1, 1.7, 3.2, 2.1, 1.7);


{ TestFFT }

{ -- compare Fast & Discrete Complex Fourier Transform -- }

procedure TTestFFT.TestComplexA(Data: PComplexArray;
                                     Count: Integer; Sign: Integer);
var
  FFTData, DFTData: TDynComplexArray;
  M: Integer;

begin
  DFTData:= ComplexDFT(Data, Count, Sign);
  FFTData:= ComplexArray(Data, Count);
  ksMath.ComplexFFT(FFTData, Count, Sign);

  for M:= 0 to Count - 1 do begin
    CheckEquals(DFTData[M].Re, FFTData[M].Re, 1E-15);
    CheckEquals(DFTData[M].Im, FFTData[M].Im, 1E-15);
  end;
end;

// "plus" complex transform test
procedure TTestFFT.TestComplexA0;
begin
  TestComplexA(@ComplexArr0, Length(ComplexArr0), 1);
  TestComplexA(@ComplexArr1, Length(ComplexArr1), 1);
  TestComplexA(@ComplexArr2, Length(ComplexArr2), 1);
end;

// "minus" complex transform test
procedure TTestFFT.TestComplexA1;
begin
  TestComplexA(@ComplexArr0, Length(ComplexArr0), -1);
  TestComplexA(@ComplexArr1, Length(ComplexArr1), -1);
  TestComplexA(@ComplexArr2, Length(ComplexArr2), -1);
end;

{ -- plus/minus Fourier Transform -- }

procedure TTestFFT.TestComplexB(Data: PComplexArray; Count: Integer);
var
  FFTData: TDynComplexArray;
  M: Integer;

begin
// plus + minus transform
  FFTData:= ComplexArray(Data, Count);
  ksMath.ComplexFFT(FFTData, Count, 0);      // plus
  ksMath.ComplexFFT(FFTData, Count, -1);     // minus
  for M:= 0 to Count - 1 do begin
    CheckEquals(Data^[M].Re, FFTData[M].Re / Count, 1E-15);
    CheckEquals(Data^[M].Im, FFTData[M].Im / Count, 1E-15);
  end;

// minus + plus transform
  FFTData:= ComplexArray(Data, Count);
  ksMath.ComplexFFT(FFTData, Count, -1);     // minus
  ksMath.ComplexFFT(FFTData, Count, 0);      // plus
  for M:= 0 to Count - 1 do begin
    CheckEquals(Data^[M].Re, FFTData[M].Re / Count, 1E-15);
    CheckEquals(Data^[M].Im, FFTData[M].Im / Count, 1E-15);
  end;
end;

procedure TTestFFT.TestComplexB0;
begin
  TestComplexB(@ComplexArr0, Length(ComplexArr0));
  TestComplexB(@ComplexArr1, Length(ComplexArr1));
  TestComplexB(@ComplexArr2, Length(ComplexArr2));
end;

{ -- real function's Fourier fransform (RealFFT) tests -- }

{ -- compare Fast & Discreet Real Fourier Transform -- }

procedure TTestFFT.TestRealA(Data: PRealArray;
                                  Count: Integer; Sign: Integer);
var
  FFTData: TDynRealArray;
  DFTData: TDynComplexArray;
  M: Integer;

begin
  Assert((Sign = -1) or (Sign = 1), 'Sign must be -1 or 1');
  DFTData:= RealDFT(Data, Count, Sign);
  FFTData:= RealArray(Data, Count);
  ksMath.RealFFT(FFTData, Count, 0);

  CheckEquals(DFTData[0].Re, FFTData[0], 1E-15);
  CheckEquals(DFTData[0].Im, 0, 1E-15);

  for M:= 1 to Count shr 1 - 1 do begin
    CheckEquals(DFTData[M].Re, FFTData[2*M], 1E-15);
    CheckEquals(DFTData[M].Im, Sign * FFTData[2*M+1], 1E-15)
  end;

  CheckEquals(DFTData[Count shr 1].Re, FFTData[1], 1E-15);
  CheckEquals(DFTData[Count shr 1].Im, 0, 1E-15);

  for M:= 1 to Count shr 1 - 1 do begin
    CheckEquals(DFTData[Count - M].Re, FFTData[2*M], 1E-15);
    CheckEquals(DFTData[Count - M].Im, - Sign * FFTData[2*M + 1], 1E-15)
  end;

end;

procedure TTestFFT.TestSinF(Func: TRealFunc; Count: Integer);
var
  DFTData, FFTData: TDynRealArray;
  N: Integer;

begin
  DFTData:= SinDFT(Func, Count);
  FFTData:= RealArray(Func, Count);
  ksMath.SinFFT(FFTData, Count);
  for N:= 0 to Count - 1 do
    CheckEquals(DFTData[N], FFTData[N], 1E-15);
  ksMath.SinFFT(FFTData, Count);
  for N:= 0 to Count - 1 do
    CheckEquals(Func(N) * Count / 2, FFTData[N], 1E-15);
end;

procedure TTestFFT.TestCos1F(Func: TRealFunc; Count: Integer);
var
  DFTData, FFTData: TDynRealArray;
  N: Integer;

begin
  DFTData:= Cos1DFT(Func, Count);
  FFTData:= RealArray(Func, Count + 1);
  ksMath.Cos1FFT(FFTData, Count);
  for N:= 0 to Count do
    CheckEquals(DFTData[N], FFTData[N], 1E-15);
  ksMath.Cos1FFT(FFTData, Count);
  for N:= 0 to Count do
    CheckEquals(Func(N) * Count / 2, FFTData[N], 1E-15);
end;

procedure TTestFFT.TestCos1A(Data: PRealArray; Count: Integer);
begin
  TestCos1F(
    function(N: Integer): Extended
    begin
      Result:= Data[N];
    end,
    Count);
end;

procedure TTestFFT.TestCos1A0;
const
  Arr1: packed array[0..4] of Extended = (1.2, 0.34, 0.55, -0.37, 0.1);
  Arr2: packed array[0..8] of Extended = (1.2, 0.34, 0.55, -0.37, 0.1,
                                         -0.4, 1.03, 0.11, -0.7);
  Arr3: packed array[0..8] of Extended = (0.2, 1.34, 0.33, -0.17, 0.3,
                                         -0.4, 1.03, 0.11, -0.7);

begin
  TestCos1A(@Arr1, Length(Arr1) - 1);
  TestCos1A(@Arr2, Length(Arr2) - 1);
  TestCos1A(@Arr3, Length(Arr3) - 1);
end;

{ === Second Form of Cosine Transform Tests === }

procedure TTestFFT.TestCos2F(Func: TRealFunc; Count: Integer);
var
  DFTData, FFTData: TDynRealArray;
  N: Integer;

begin
                                      // test forward transform
  DFTData:= Cos2DFT(Func, Count, 0);
  FFTData:= RealArray(Func, Count);
  ksMath.Cos2FFT(FFTData, Count, 0);
  for N:= 0 to Count - 1 do
    CheckEquals(DFTData[N], FFTData[N], 1E-15);
                                      // test inverse transform
  ksMath.Cos2FFT(FFTData, Count, -1);
  for N:= 0 to Count - 1 do
    CheckEquals(Func(N) * Count / 2, FFTData[N], 1E-15);
end;

procedure TTestFFT.TestCos2A(Data: PRealArray; Count: Integer);
begin
  TestCos2F(
    function(N: Integer): Extended
    begin
      Result:= Data[N];
    end,
    Count);
end;

procedure TTestFFT.TestCos2A0;
const
  Arr1: packed array[0..3] of Extended = (1.2, 0.34, 0.55, -0.37);
  Arr2: packed array[0..7] of Extended = (1.2, 0.34, 0.55, -0.37, 0.1,
                                         -0.4, 1.03, 0.11);
  Arr3: packed array[0..7] of Extended = (0.2, 1.34, 0.33, -0.17, 0.3,
                                         -0.4, 1.03, 0.11);

begin
  TestCos2A(@Arr1, Length(Arr1));
  TestCos2A(@Arr2, Length(Arr2));
  TestCos2A(@Arr3, Length(Arr3));
end;

procedure TTestFFT.TestSinF0;
begin
  TestSinF( function(N: Integer): Extended
    begin
      Result:= Sin(pi / 4 * N) + 0.5 * Sin(pi / 2 * N);
    end,
    8);
  TestSinF( function(N: Integer): Extended
    begin
      Result:= Sin(pi / 8 * N) + 0.25 * Sin(pi / 2 * N);
    end,
    16);
  TestSinF( function(N: Integer): Extended
    begin
      Result:= Sin(pi / 8 * N) - 0.5 * Sin(pi / 4 * N);
    end,
    16);
  TestSinF( function(N: Integer): Extended
    begin
      Result:= Sin(pi / 32 * N) - 0.2 * Sin(pi / 24 * N);
    end,
    8);
end;

procedure TTestFFT.TestRealB(Data: PRealArray; Count, Sign: Integer);
var
  FFTData: TDynRealArray;
  DFTData: TDynComplexArray;
  M: Integer;

begin
  Assert((Sign = -1) or (Sign = 1), 'Sign must be -1 or 1');
  DFTData:= RealDFT(Data, Count, Sign);

  SetLength(FFTData, Count);
  FFTData[0]:= DFTData[0].Re;
  FFTData[1]:= DFTData[Count shr 1].Re;

  for M:= 1 to (Count shr 1) - 1 do begin
    FFTData[2*M]:= DFTData[M].Re;
    FFTData[2*M+1]:= Sign * DFTData[M].Im;
  end;

  ksMath.RealFFT(FFTData, Count, -1);

  for M:= 0 to Count - 1 do
    CheckEquals(Data[M] * (Count shr 1), FFTData[M], 1E-15);
end;

procedure TTestFFT.TestRealC(Data: PRealArray; Count: Integer);
var
  FFTData: TDynRealArray;
  M: Integer;

begin
  FFTData:= RealArray(Data, Count);
  ksMath.RealFFT(FFTData, Count, 0);
  ksMath.RealFFT(FFTData, Count, -1);
  for M:= 0 to Count - 1 do
    CheckEquals(Data[M] * (Count shr 1), FFTData[M], 1E-15);

  FFTData:= RealArray(Data, Count);
  ksMath.RealFFT(FFTData, Count, -1);
  ksMath.RealFFT(FFTData, Count, 0);
  for M:= 0 to Count - 1 do
    CheckEquals(Data[M] * (Count shr 1), FFTData[M], 1E-15);
end;

procedure TTestFFT.TestRealF(Func: TRealFunc; Count: Integer; Sign: Integer);
var
  FFTData: TDynRealArray;
  DFTData: TDynComplexArray;
  N: Integer;

begin
  Assert((Sign = -1) or (Sign = 1), 'Sign must be -1 or 1');
  FFTData:= RealArray(Func, Count);
  DFTData:= RealDFT(@FFTData[0], Count, Sign);
  ksMath.RealFFT(@FFTData[0], Count, 0);

  CheckEquals(FFTData[0], DFTData[0].Re, 1E-15);
  CheckEquals(FFTData[1], DFTData[Count shr 1 - 1].Re, 1E-15);
  CheckEquals(0, DFTData[0].Im, 1E-15);
  CheckEquals(0, DFTData[Count shr 1 - 1].Im, 1E-15);

  for N:= 1 to Count shr 1 - 1 do begin
    CheckEquals(DFTData[N].Re, FFTData[2 * N], 1E-15);
    CheckEquals(Sign * DFTData[N].Im, FFTData[2 * N + 1], 1E-15);
    CheckEquals(DFTData[Count - N].Re, FFTData[2 * N], 1E-15);
    CheckEquals(-Sign* DFTData[Count - N].Im, FFTData[2 * N + 1], 1E-15);
  end;
end;

procedure TTestFFT.TestRealA0;
begin
  TestRealA(@RealArr0, Length(RealArr0), 1);
  TestRealA(@RealArr1, Length(RealArr1), 1);
  TestRealA(@RealArr2, Length(RealArr2), 1);
end;

procedure TTestFFT.TestRealA1;
begin
  TestRealA(@RealArr0, Length(RealArr0), -1);
  TestRealA(@RealArr1, Length(RealArr1), -1);
  TestRealA(@RealArr2, Length(RealArr2), -1);
end;

procedure TTestFFT.TestRealB0;
begin
  TestRealB(@RealArr0, Length(RealArr0), 1);
  TestRealB(@RealArr1, Length(RealArr1), 1);
  TestRealB(@RealArr2, Length(RealArr2), 1);
end;

procedure TTestFFT.TestRealB1;
begin
  TestRealB(@RealArr0, Length(RealArr0), -1);
  TestRealB(@RealArr1, Length(RealArr1), -1);
  TestRealB(@RealArr2, Length(RealArr2), -1);
end;

procedure TTestFFT.TestRealC0;
begin
  TestRealC(@RealArr0, Length(RealArr0));
  TestRealC(@RealArr1, Length(RealArr1));
  TestRealC(@RealArr2, Length(RealArr2));
end;

procedure TTestFFT.TestRealF0;
begin
  TestRealF(FFunc0, 8, 1);
  TestRealF(FFunc1, 16, 1);
  TestRealF(FFunc2, 16, 1);
  TestRealF(FFunc0, 8, -1);
  TestRealF(FFunc1, 16, -1);
  TestRealF(FFunc2, 16, -1);
end;

procedure TTestFFT.TestSpectrumA(Data: PRealArray; Count: Integer);
begin
  TestSpectrumF(
    function(N: Integer): Extended
    begin
      Result:= Data[N];
    end,
    Count
  );
end;

procedure TTestFFT.TestSpectrumA0;
begin
  TestSpectrumA(@RealArr0, Length(RealArr0));
  TestSpectrumA(@RealArr1, Length(RealArr1));
  TestSpectrumA(@RealArr2, Length(RealArr2));
end;

procedure TTestFFT.TestSpectrumF(Func: TRealFunc; Count: Integer);
var
  FFTData: TDynRealArray;
  DFTData: TDynComplexArray;
  N: Integer;
  Value: TksComplex;

begin
  DFTData:= RealDFT(Func, Count, -1);

  FFTData:= RealArray(Func, Count);
  ksMath.RealFFT(FFTData, Count, 0);

  for N:= 0 to Count - 1 do begin
    Value:= RFTSpectrum(FFTData, Count, N);
    CheckEquals(DFTData[N].Re, Value.Re, 1E-15);
    CheckEquals(DFTData[N].Im, Value.Im, 1E-15);
  end;
end;

procedure TTestFFT.TestSpectrumF0;
begin
  TestSpectrumF(FFunc0, 8);
  TestSpectrumF(FFunc1, 16);
  TestSpectrumF(FFunc2, 16);
end;

procedure TTestFFT.TestTwoA(Data1, Data2: PRealArray; Count: Integer);
var
  DFTData1, DFTData2: TDynComplexArray;
  FFTData1, FFTData2: TDynComplexArray;
  N: Integer;

begin
  DFTData1:= RealDFT(Data1, Count, 0);
  DFTData2:= RealDFT(Data2, Count, 0);
  SetLength(FFTData1, Count);
  SetLength(FFTData2, Count);
  ksMath.RealFFT2(Data1, Data2, FFTData1, FFTData2, Count);
  for N:= 0 to Count - 1 do begin
    CheckEquals(DFTData1[N].Re, FFTData1[N].Re, 1E-15);
    CheckEquals(DFTData1[N].Im, FFTData1[N].Im, 1E-15);

    CheckEquals(DFTData2[N].Re, FFTData2[N].Re, 1E-15);
    CheckEquals(DFTData2[N].Im, FFTData2[N].Im, 1E-15);
  end;
end;

procedure TTestFFT.TestTwoA0;
begin
  TestTwoA(@RealArr1, @RealArr2, 8);
end;

procedure TTestFFT.TestTwoF(Func1, Func2: TRealFunc; Count: Integer);
var
  Data1, Data2: TDynRealArray;
  DFTData1, DFTData2: TDynComplexArray;
  FFTData1, FFTData2: TDynComplexArray;
  N: Integer;

begin
  Data1:= RealArray(Func1, Count);
  Data2:= RealArray(Func2, Count);

  DFTData1:= RealDFT(Func1, Count, 0);
  DFTData2:= RealDFT(Func2, Count, 0);

  SetLength(FFTData1, Count);
  SetLength(FFTData2, Count);

// Two - function fast Fourier transform
  RealFFT2(@Data1[0], @Data2[0], @FFTData1[0], @FFTData2[0], Count);

  for N:= 0 to Count - 1 do begin
    CheckEquals(DFTData1[N].Re, FFTData1[N].Re, 1E-15);
    CheckEquals(DFTData1[N].Im, FFTData1[N].Im, 1E-15);

    CheckEquals(DFTData2[N].Re, FFTData2[N].Re, 1E-15);
    CheckEquals(DFTData2[N].Im, FFTData2[N].Im, 1E-15);
  end;
end;

procedure TTestFFT.TestTwoF0;
begin
  TestTwoF(FFunc1, FFunc2, 16);
end;

procedure TTestFFT.TestTwoG(Func1, Func2: TRealFunc; Count: Integer);
var
  Data1, Data2: TDynRealArray;
  FFTData1, FFTData2: TDynComplexArray;
  N: Integer;

begin
  Data1:= RealArray(Func1, Count);
  Data2:= RealArray(Func2, Count);

  SetLength(FFTData1, Count);
  SetLength(FFTData2, Count);

// Two - function fast Fourier transform
  RealFFT2(@Data1[0], @Data2[0], @FFTData1[0], @FFTData2[0], Count);

// fast Fourier transform
  RealFFT(@Data1[0], Count, 0);
  RealFFT(@Data2[0], Count, 0);

  CheckEquals(Data1[0], FFTData1[0].Re, 1E-15);
  CheckEquals(Data1[1], FFTData1[Count shr 1 - 1].Re, 1E-15);
  CheckEquals(0, FFTData1[0].Im, 1E-15);
  CheckEquals(0, FFTData1[Count shr 1 - 1].Im, 1E-15);

  CheckEquals(Data2[0], FFTData2[0].Re, 1E-15);
  CheckEquals(Data2[1], FFTData2[Count shr 1 - 1].Re, 1E-15);
  CheckEquals(0, FFTData2[0].Im, 1E-15);
  CheckEquals(0, FFTData2[Count shr 1 - 1].Im, 1E-15);

  for N:= 1 to Count shr 1 - 1 do begin
    CheckEquals(Data1[2 * N], FFTData1[N].Re, 1E-15);
    CheckEquals(Data1[2 * N + 1], FFTData1[N].Im, 1E-15);
    CheckEquals(Data1[2 * N], FFTData1[Count - N].Re, 1E-15);
    CheckEquals(Data1[2 * N + 1], - FFTData1[Count - N].Im, 1E-15);

    CheckEquals(Data2[2 * N], FFTData2[N].Re, 1E-15);
    CheckEquals(Data2[2 * N + 1], FFTData2[N].Im, 1E-15);
    CheckEquals(Data2[2 * N], FFTData2[Count - N].Re, 1E-15);
    CheckEquals(Data2[2 * N + 1], - FFTData2[Count - N].Im, 1E-15);
  end;
end;

procedure TTestFFT.TestTwoG0;
begin
  TestTwoG(FFunc1, FFunc2, 16);
end;

procedure TTestFFT.TestCorrA(Data1, Data2: PRealArray; Count: Integer);
var
  Corr1, Corr2: TDynRealArray;
  N: Integer;

begin
  Corr1:= DirectRealCorr(Data1, Data2, Count);
  SetLength(Corr2, Count);
  RealCorr(Data1, Data2, Corr2, Count);
  for N:= 0 to Count - 1 do
    CheckEquals(Corr1[N], Corr2[N], 1E-15);
end;

procedure TTestFFT.TestCorrA0;
begin
  TestCorrA(@RealArr1, @RealArr2, 8);
end;

procedure TTestFFT.TestCorrF(Func1, Func2: TRealFunc; Count: Integer);
var
  Data1, Data2: TDynRealArray;
  Corr1, Corr2: TDynRealArray;
  N: Integer;

begin
  Corr1:= DirectRealCorr(Func1, Func2, Count);

  Data1:= RealArray(Func1, Count);
  Data2:= RealArray(Func2, Count);
  SetLength(Corr2, Count);
  RealCorr(Data1, Data2, Corr2, Count);

  for N:= 0 to Count - 1 do
    CheckEquals(Corr1[N], Corr2[N], 1E-15);
end;

procedure TTestFFT.TestCorrF0;
begin
  TestCorrF(FFunc1, FFunc2, 16);
end;

{ -- real function's circular autocorrelation (RealAutoCorr) tests -- }

procedure TTestFFT.TestAutoCorrA(Data: PRealArray; Count: Integer);
begin
  TestAutoCorrF(
    function(N: Integer): Extended
    begin
      Result:= Data[N];
    end,
    Count
  );
end;

procedure TTestFFT.TestAutoCorrA0;
begin
  TestAutoCorrA(@RealArr0, Length(RealArr0));
  TestAutoCorrA(@RealArr1, Length(RealArr1));
  TestAutoCorrA(@RealArr2, Length(RealArr2));
end;

procedure TTestFFT.TestAutoCorrF(Func: TRealFunc; Count: Integer);
var
  Corr1, Corr2: TDynRealArray;
  N: Integer;

begin
  Corr1:= DirectRealAutoCorr(Func, Count);
  Corr2:= RealArray(Func, Count);
  RealAutoCorr(Corr2, Count);

  for N:= 0 to Count - 1 do
    CheckEquals(Corr1[N], Corr2[N], 1E-15);
end;

procedure TTestFFT.TestAutoCorrF0;
begin
  TestAutoCorrF(FFunc0, 8);
  TestAutoCorrF(FFunc1, 16);
  TestAutoCorrF(FFunc2, 16);
end;

procedure TTestFFT.SetUp;
begin
  FFunc0:= function(N: Integer): Extended
    begin
      Result:= Sin(pi / 4 * N) + 0.5 * Sin(pi / 2 * N + 3 * pi / 4);
    end;

  FFunc1:= function(N: Integer): Extended
    begin
      Result:= Sin(pi / 8 * N) + 0.25 * Sin(pi / 2 * N + 3 * pi / 8);
    end;

  FFunc2:= function(N: Integer): Extended
    begin
      Result:= Sin(pi / 8 * N) - 0.5 * Sin(pi / 4 * N + pi / 8);
    end;
end;

initialization
  // Register any test cases with the test runner
  RegisterTest(TTestFFT.Suite);
end.

